#include <complex>
#include <iostream>
#include <tnt.h>
typedef std::complex<double> DComplex;
typedef TNT::Array2D<DComplex> ComplexMatrix;
typedef TNT::Array1D<DComplex> ComplexArray;
typedef TNT::Array2D<double > DoubleMatrix;
typedef TNT::Array1D<double > DoubleArray;


typedef int (* SELECT)(std::complex<double> *);

extern "C" void 
zgees_(char *jobvs, char *sort, SELECT select, int *n,
       DComplex *a, int *lda, int *sdim, DComplex *w, DComplex *vs, int *ldvs,
       DComplex *work, int *lwork, double *rwork, int *bwork, int *info);

void 
zgees(char jobvs, char sort, SELECT select, int n, DComplex *a, int lda,
      int &sdim, DComplex *w, DComplex *vs, int ldvs, DComplex *work,
      int lwork, double *rwork, int *bwork, int &info)
{
  
  zgees_(&jobvs, &sort, select, &n, a, &lda, &sdim, w, vs, &ldvs, work, &lwork,
	rwork, bwork, &info);
}




double amatrix[] = {-2.13, 0.51, 0.69,-0.04,-1.05,
		    0.81,-1.41,-0.93,-1.40,-0.25,
		    1.05,-0.59,-0.66,-1.50,-0.04,
		    -1.02, 0.39, 1.32,-1.43, 1.40,
		    0.69,-0.52,-0.04,-0.52, 0.68};


int main(){
  const int N = 5;
  DoubleMatrix A(N, N, amatrix);
  ComplexMatrix Ac(N, N);
  
  for (int i = 0; i < N; ++i)
    for (int j = 0; j < N; j++)
      Ac[i][j] = std::complex<double>(A[i][j], 0.0);
  int sdim, info;
  ComplexArray W(N);
  int lwork = 4*N;
  ComplexArray work(lwork);
  DoubleArray rwork(N);
  
  zgees('N', 'N', 0, N, Ac[0], N, sdim, W, 0, 1, work, lwork, rwork, 0, info);
  using namespace std;
  
  cout << Ac;
  
  
}


